# Dependency: sim.results in "sim_results_miss1.Rdata"

library(reshape)
library(ggplot2)

# Load simulation output

load("sim_results_miss1.Rdata")

J <- 1000
K<- 17

# Bias formula

bias <- function(v,true) {
  sum(v - true)/length(v)
}

true.effect <- NULL
sim.results.long <- NULL
sim.results.mtx <- NULL
est.bias <- NULL

for (k in 1:K) {

true.effect[[k]] <- sim.results.miss[[k]][1]
sim.results.long[[k]] <- sim.results.miss[[k]][-1]

sim.results.mtx[[k]] <- matrix(NA, nrow = length(sim.results.long[[k]]), ncol = J)

for (i in 1:length(sim.results.long[[k]])) { # ITERATE OVER SIMULATIONS
  
  sim.results.mtx[[k]][i,] <- as.matrix(sim.results.long[[k]][[i]][,1])
  
}

est.bias[[k]] <- apply(sim.results.mtx[[k]], 1, bias, true = true.effect[[k]][[1]][[2]])

names(est.bias[[k]]) <- names(sim.results.long[[k]])

}

# Plot bias by value of heterogeneous effect of G

est.bias.matrix <- matrix(unlist(est.bias), byrow = T, ncol = length(sim.results.long[[1]]))

colnames(est.bias.matrix) <- names(sim.results.long[[1]])

betahg.vector <- c(-1.6, -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6) 

rownames(est.bias.matrix) <- betahg.vector

weighted.methods <- c("vATT.wweight", "vATT.wweight.reg", "vATT.wnn", "vATT.wnn.reg", "vATT.ws", "vATT.ws.reg")

data.plot <- melt(est.bias.matrix[,colnames(est.bias.matrix) %in% weighted.methods])

names(data.plot) <- c("Beta_G", "method", "bias")

# Figure 2, Panel A

ggplot(data=data.plot, aes(x=Beta_G, y=abs(bias), group=method)) +
  geom_line() + labs( x = expression(Beta[S[2]]))  + expand_limits(y=c(0,0.4)) + theme_bw(base_size = 18)